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i We present a model-independent investigation of the Wilkinson Microwave 

Anisotropy Probe (WMAP) data with respect to scale-independent and scale- 

£j dependent non-Gaussianities (NGs). To this end, we employ the method of con- 

strained randomization. For generating so-called surrogate maps a well-specified shuf- 

i i fling scheme is applied to the Fourier phases of the original data, which allows to test 

for the presence of higher order correlations (HOCs) also and especially on well-defined 
scales. 

Using scaling indices as test statistics for the HOCs in the maps we find highly 
significant signatures for non-Gaussianities when considering all scales. We test for 
NGs in four different /-bands A/, namely in the bands Al = [2,20], Al = [20,60], 
^ Al = [60,120] and Al = [120,300]. We find highly significant signatures for both 

^ . non-Gaussianities and ecliptic hemispherical asymmetries for the interval Al = [2, 20] 

covering the large scales. We also obtain highly significant deviations from Gaussian- 
ity for the band Al = [120,300]. The result for the full /-range can then easily be 
interpreted as a superposition of the signatures found in the bands Al = [2,20] and 
Al = [120,300]. We find remarkably similar results when analyzing different ILC-like 
^ maps based on the WMAP three, five and seven year data. We perform a set of tests 

to investigate whether and to what extend the detected anomalies can be explained by 
systematics. While none of these tests can convincingly rule out the intrinsic nature 
of the anomalies for the low I case, the ILC map making procedure and/or residual 
noise in the maps can also lead to NGs at small scales. 

Our investigations prove that there are phase correlations in the WMAP data of the 
CMB. In the absence of an explanation in terms of Galactic foregrounds or known 
systematic artefacts, the signatures at low I must so far be taken to be cosmological at 
high significance. These findings would strongly disagree with predictions of isotropic 
cosmologies with single field slow roll inflation. 

The task is now to elucidate the origin of the phase correlations and to understand 
the physical processes leading to these scale-dependent non-Gaussianities - if it turns 
out that systematics as cause for them must be ruled out. 
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The Cosmic Microwave Background (CMB) radiation rep- 
resents the oldest observable signal in the Universe. Since 
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this relic radiation has its origin just 380000 years after 
the Big Bang when the CMB photons were last scattered 
off electrons, this radiation is one of the most important 
sources of information to gain more knowledge about the 
very early Universe. Estimating the linear correlations of 
the temperature fluctuations in the CMB as measured e.g. 
with the WMAP satellite by means of the power spectrum 
has yielded very precise determinations of the parameters 
of the standard ACDM cosmological model like the age, the 
geometry and the matter and energy content of the Universe 
( Komatsu et al. ||20Q9a"l |2011 ). 

Analyzing CMB maps by means of the power spectrum rep- 
resents an enormous compression of information contained 
in the data from approx. 10 6 temperature values to roughly 
1000 numbers for the power spectrum. It has often been 
pointed out ( |Komatsu et al. ( 2009b| and references therein) 
that this data compression is lossless and thus fully jus- 
tified, if and only if the statistical distribution of the ob- 
served fluctuations is a Gaussian distribution with random 
phases. Any information that is contained in the phases and 
the correlations among them, is not encoded in the power 
spectrum, but has to be extracted from measurements of 
higher-order correlation (HOC). Thus, the presence of phase 
correlations may be considered as an unambiguous evidence 
of non-Gaussianity (NG). Otherwise, non-Gaussianity can 
only be defined by the negation of Gaussianity. 
Primordial NG represents one way to test theories of in- 
flation with the ultimate goal to constrain the shape of 
the potential of the inflaton field (s) and their possible 
(self-) interactions. While the simplest single field slow roll 
inflationary scenario predicts that fluctuations are nearly 
Gaussian (|Guth||1981| iLindel 119821 lAlbrecht & Steinhardl 



et al.| fl2004|); |Hansen, Banday, fe G6rski| pXte); |Vielva et 



1982]), a variety of more complex models predict deviations 



from Gaussianity (Linde & Mukhanov ||1997 


Peebles 1997 


Bernardeau & Uzan 2002 Acquaviva et al. 


2003). Models 



in which the Lagrangian is a general function of the infla- 



ton and powers of its first derivative ( Armendariz-Picon, 



Damour fc Mukhanov 1999 Garriga & Mukhanov 1999) 



can lead to scale-dependent non-Gaussianities, if the sound 
speed varies during inflation. Similarly, string theory mod- 
els that give rise to large non-Gaussianity have a natural 
scale dependence ( |Chen |2005||Lo Verde et al. |2008| . Also, 
NGs put strong constraints on alternatives to the inflation- 
ary paradigm (Buchbinder7^hour^^Ovrur]20^ Lehners 



& Steinhardt 12008). 



Given the plethora of conceivable scenarios for the very early 
Universe, it is worth first checking what is in the data in a 
model-independent way. Further, such a model-independent 
approach has a large discovery potential to detect yet un- 
expected fingerprints of nonlinear physics in the early uni- 
verse. Thus, a detection of possibly scale-dependent non- 
Gaussianity being encrypted in the phase correlations in the 
WMAP data would be of great interest. While a detection of 
non-Gaussianity could be indicative of an experimental sys- 
tematic effect or of residual foregrounds, it could also point 
to new cosmological physics. 

The investigations of deviations from Gaussianity in the 
CMB (see |Komatsu et al. | ( |2009a| and references therein) 
and claims for the detection of non-Gaussianitiy and a vari- 
ety of other anomalies like hemispherical asymmetries, lack 
of power at large angular scales, alignment of multipoles, de- 



(2004) 



Eriksen et al. (20051 2007); de Oliveira- Costa et 



2004); Rath et al. (2007); McEwen et al. (2008); Ross- 



manith et al. (2009); Copi et al. (2009 2010); Hansen et al. 



( 2009| ); |Yoho, Ferrer fe Starkman | jf2010)) have been made, 
where the statistical significance of some of the detected 



signatures is still subject to discussion (Zhang & Huterer 
|2010| |Bennett et al.||2011| ) . These studies have in common 
that the level of non-Gaussianity is assessed by comparing 
the results for the measured data with simulated CMB- 
maps which were generated on the basis of the standard 
cosmological model and/or specific assumptions about the 
nature of the non-Gaussianities as parametrized with e.g. 
the scalar, scale-independent parameter f n i. Other studies 
focused on the detection of signatures in the distribution of 
Fourier phases flChiang et al.|2003||Coles et al.|2004| |Nasel- 
|sky e t al. 2005j |Chiang, Naselsky &; Coles |2007| representing 
deviations from the random phase hypothesis for Gaussian 
random fields. These model-independent tests also revealed 
signatures of NGs. Pursuing this approach one can go one 
step further and investigate possible phase correlations and 
their relation to the morphology of the CMB maps by means 
of so-called surrogate maps. 

This technique of surrogate data sets (Theil er"et al. | |l992) 
was originally developed for nonlinear time series analysis. 
In this field of research complex systems like the climate, 
stock-market, heart-beat variability, etc. are analyzed (see 



e.g. |Bunde, Kropp &; Schellnhuber ( 2002| ) and references 
therein). For those systems a full modeling is barely or not 
possible. Therefore, statistical methods of constrained ran- 
domization involving surrogate data sets were developed to 
infer some information about the nature of the underly- 
ing physical process in a completely data-driven, i.e. model- 
independent way. One of the first and most basic question 
here is whether a (quasiperiodic) process is completely lin- 
ear or whether also weak nonlinearities can be detected in 
the data. The basic formalism to answer this question is to 
compute statistics sensitive to HOCs for the original data 
set and for an ensemble of surrogate data sets, which mimic 
the linear properties of the original time series while wiping 
out all phase correlations. If the computed measure for the 
original data is significantly different from the values ob- 
tained for the set of surrogates, one can infer that the data 
contain HOCs. 

Extensions of this formalism to three-dimensional galaxy 
distributions (Rat h et al.||2002| and two -dimensional sim- 
ulated flat CMB maps ( Rath &; Schuecker||2003 ) have been 
proposed and discussed. By introducing a more sophisti- 
cated two-step surrogatization scheme for full-sky CMB ob- 
servations it has become possible to also test for scale- 



dependent NG in a model-independent way ( |Rath et al. 



2009 ). Probing NG on the largest scales (I < 20) yielded 
highly significant signatures for both NG and ecliptic hemi- 
spherical asymmetries. 

In this paper, we apply the method of constrained random- 
ization to the WMAP five year and seven year data in or- 
der to test for scale-independent and scale-dependent non- 
Gaussianity up to / = 300 as encoded in the Fourier phase 
correlations. Further, this work fully recognises the need to 
rule out foregrounds and systematic artefacts as the ori- 



gin of the detections (as advised by Bennett et al. (2011)) 



tection of the Cold Spot etc. (see e.g. Park (2004); Eriksen Therefore, a large part of our analyses is dedicated to var- 
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ious checks on systematics to single out possible causes of 
the detected anomalies. 

The paper is organized as follows: In Section 2 we briefly 
describe the observational and simulated data we use in our 
study. The method of constrained randomization is reviewed 
in some detail in Section 3. Scaling indices, which we use as 
test statistic, and the statistics derived out of them are dis- 
cussed in Section 4. In Section 5 we present our results and 
we draw our conclusions in Section 6. 



2 DATA SETS 

We used the seven years foreground-cleaned internal linear 



combination (ILC) map (Gold_et al. 2011) generated and 



provided by the WMAP tearrQ (in the following: ILC7) 
For comparison we also included the map produced by |De- 



labrouille et al. (2009), namely the five years needlet-based 



ILC map, which has been shown to be significantly less con- 
taminated by foreground and noise than other existing maps 
obtained from WMAP data (in the following: NILC5). 
To check for systematics we also analyzed the following set 
of maps: 

1) Uncorrected ILC map 

The ILC map is a weighted linear combination of the 5 fre- 
quency channels that recovers the CMB signal. The weights 
are derived by requiring minimum variance in a given re- 
gion of the sky under the constraint that the sum of the 
weights is unity. Such weights, however, cannot null an ar- 
bitrary foreground signal with a non-blackbody frequency 
spectrum, thus some residuals due to Galactic emission will 
remain. The WMAP team attempts to correct for this "bias" 
with an estimation of the residual signal based on simula- 
tions and a model of the foreground sky. Our uncorrected 
map (UILC7 in the following) is simply the ILC without ap- 
plying this correction, computed from the weights provided 
|Gold et al.~l fl2011| ) and the 1-degree smoothed WMAP 



data. 

2) Asymmetric beam map 

Beam asymmetries may result in statistically anisotropic 
CMB maps . To asses these effects on the signatures of scale- 
dependent NGs and their (an-)isotropies we make use of the 
publicly available CMB sky simulations including the effects 
of asymmetric beams (Wehus et al. 2009). Specifically we 
analyse a simulated map of the Vl-band, because this band 
is considered to have the least foreground contamination. 

3) Simulated coadded VW-band map 

To make sure that neither systematic effects are induced by 
the method of constrained randomization nor the WMAP- 
like beam and noise properties lead to systematic deviations 
from Gaussianity we include in our analysis a co-added VW- 
map as obtained using the standard ACDM best fit power 
spectrum and WMAP-like beam and noise properties. Note 
that this map did not undergo the ILC-map making proce- 
dure. 

4) Simulated ILC map 

Simulated sky maps result from processing a simulated dif- 
ferential time-ordered data (TOD) stream through the same 
calibration and analysis pipeline that is used for the flight 



http:/ /lambda, gsfc.nasa.gov 



data. The TOD is generated by sampling a reference sky 
that includes both CMB and Galactic foregrounds with 
the actual flight pointings, and adding various instrumen- 
tal artefacts. We have then processed the individual result- 
ing data into 7 separate simulated yearly ILC maps, plus a 
7-year merge. It is worth noting that, if the yearly frequency- 
averaged maps are combined into ILCs using the |Gold etaL] 
H fl20TT| ) 7 -year weights per region, then the resulting ILCs 
show clear Galactic plane residuals. This reflects the fact 
that the simulated data has a different CMB realisation to 
the observed sky, and may additionally represent a mismatch 
between the simulated foreground properties and the true 
sky in the Galactic plane. Instead, we analyse the 7- year 
merged simulated data to compute the ILC weights for the 
simulations, then apply to all yearly data sets separately. 
However, the derived weights are quite different from the 
WMAP7 ones, which would imply different noise properties 
in the simulated ILC data compared to the real data. Care 
should be exercised for any results that are sensitive to the 
specific noise pattern. 
5) Difference ILC map 

Finally, we consider the difference map (year 7 - year 6) from 
yearly ILC-maps computed using the same weights and re- 



gions as the 7- year data set from |Gold et al. | ( |201l| ). No de- 
biasing has been applied. With this map we estimate what 
effect possible ILC-residuals may have on the detection of 
NGs. 



3 GENERATING SURROGATE MAPS 

To test for scale-dependent non-Gaussianities in a model- 
independent way we apply a two-step procedure that has 

120091). Let us 



been proposed and discussed inlRath et al. 



describe the various steps for generating surrogate maps in 
more detail: 

Consider a CMB map T(#,0), where T(#, 0) is Gaussian 
distributed and its Fourier transform. The Fourier coeffi- 
cients aim can be written as ai m — \ai m \e tcf)lrn with <^ m = 
arctan (Im(a^ rn )/Re(a^ rn )). The linear or Gaussian proper- 
ties of the underlying random field are contained in the abso- 
lute values \ai m | , whereas all HOCs - if present - are encoded 
in the phases <^i m and the correlations among them. Having 
this in mind, a versatile approach for testing for scale depen- 
dent non-Gaussianities relies on a scale-dependent shuffling 
procedure of the phase correlations followed by a statistical 
comparison of the so-generated surrogate maps. 
However, the Gaussianity of the temperature distribution 
and the randomness of the set of Fourier phases in the 
sense that they are uniformly distributed in the interval 
[— 7r,7r], are a necessary prerequisite for the application of 
the surrogate- generating algorithm, which we propose in the 
following. To fulfill these two conditions, we perform the fol- 
lowing preprocessing steps. First, the maps are remapped 
onto a Gaussian distribution in a rank-ordered way. This 
means that the amplitude distribution of the original tem- 
perature map in real space is replaced by a Gaussian distri- 
bution in a way that the rank-ordering is preserved, i.e. the 
lowest value of the original distribution is replaced with the 
lowest value of the Gaussian distribution etc. By applying 
this remapping we automatically focus on HOCs induced by 
the spatial correlations in the data while excluding any ef- 
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Figure 1. Deviations S((a(rk))) of the rotated hemispheres for three scales r^, k = 2, 6, 10 (from top to bottom) for the ILC7 map and 
for (from left to right) the shuffling intervals Al = [2, 1024], Al = [2, 20], Al = [20, 60], Al = [60, 120] and Al = [120, 300]. The expected 
correspondence between the shuffling range Al and the scales of the scale-dependent higher order statistics (a(rfc)), for which the 
largest deviations are detected, becomes apparent. While the ecliptic hemispherical asymmetries for Al = [2, 20] are most pronounced 
for the largest scaling range no (second column), the deviation S becomes largest for n when shuffling the phases of the smallest scales 
Al = [120,300] (rightmost column). 



fects coming from deviations of the temperature distribution 
from a Gaussian one. 

To ensure the randomness of the set of Fourier phases we 
performed a rank-ordered remapping of the phases onto a set 
of uniformly distributed ones followed by an inverse Fourier 
transformation. These two preprocessing steps only have 
marginal influence to the maps. The main effect is that the 
outliers in the temperature distribution are removed. Due 
to the large number of temperature values (and phases) we 
did not find any significant dependence of the specific Gaus- 
sian (uniform) realization used for remapping of the tem- 
peratures (phases) . The resulting maps may already be con- 
sidered as a surrogate map and we named it zeroth order 
surrogate map. The first and second order surrogate maps 
are obtained as follows: 

We first generate a first order surrogate map, in which any 
phase correlations for the scales, which are not of interest, 
are randomized. This is achieved by a random shuffle of the 
phases 0j m for I £ Al = [lmin, Imax], < m < / and by per- 
forming an inverse Fourier transformation. 
In a second step, N (N = 500 throughout this study) re- 
alizations of second order surrogate maps are generated for 
the first order surrogate map, in which the remaining phases 
(j)im with / G A/,0 < ra ^ / are shuffled, while the already 
randomized phases for the scales, which are not under con- 
sideration, are preserved. Note that the Gaussian properties 
of the maps, which are given by \ai m \, are exactly preserved 
in all surrogate maps. 

So far, we have applied the method of surrogates only 
to the /-range Al = [2,20]. In this paper we will repeat 
the investigations for this /-interval but using newer CMB 
maps. Furthermore, we extend the analysis to smaller scales. 
Namely, we consider three more /-intervals Al — [20,60], 
Al = [60, 120] and Al = [120, 300]. The choice of 60 as l min 
and Imax is somewhat arbitrary, whereas the lmin — 120 
and Imax = 300 for the last /-interval was selected in such 
a way that the first peak in the power spectrum is covered. 
Going to even higher /'s doesn't make much sense, because 
the ILC7 map is smoothed to 1 degree FWHM. Some other 



maps which we included in our study - especially NILC5 
- are not smoothed and we could in principle go to higher 
/'s. But to allow for a consistent comparison of the results 
obtained with the different observed and simulated input 
maps we restrict ourselves to only investigate /-intervals up 
to Imax = 300 in this study. 

Besides this two-step procedure aiming at a dedicated scale- 
dependent search of non-Gaussianity, we also test for non- 
Gaussianity using surrogate maps without specifying certain 
scales. In this case there are no scales, which are not of inter- 
est, and the first step in the surrogate map making procedure 
becomes dispensable. The zeroth order surrogate map is to 
be considered here as first order surrogate and the second 
order surrogates are generated by shuffling all phases with 
< m ^ / for all available Z's, i.e. in our case Al = [2, 1024]. 
Finally, for calculating scaling indices to test for higher 
order correlations the surrogate maps were degraded to 
N S ide = 256 and residual monopole and dipole contribu- 
tions were subtracted. The statistical comparison of the two 
classes of surrogates will reveal, whether possible HOCs on 
certain scales have left traces in the first order surrogate 
maps, which were then deleted in the second order surro- 
gates. Before the results of such a comparison of the surro- 
gate maps are shown in detail, we review the formalism of 
scaling indices. 



4 WEIGHTED SCALING INDICES AND TEST 
STATISTICS 

As test statistics for detecting and assessing possible scale- 
dependent non-Gaussianities in the CMB data weighted 



scaling indices are calculated (Rath et al. 2002 Rath &; 



Schuecker 2003 ) . The basic ideas of the scaling index method 



(SIM) stem from the calculation of the dimensions of at- 
tractors in nonlinear time series analysis (Grassberger & 



|Procacci~ 1983). Scaling indices essentially represent one 



way to estimate the local scaling properties of a point set in 
an arbitrary d-dimensional embedding space. The technique 
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offers the possibility of revealing local structural character- 
istics of a given point distribution. Thus, point-like, string- 
like and sheet-like structures can be discriminated from each 
other and from a random background. The alignment of e.g. 
string-like structures can be detected by using a proper met- 



ric for calculating the distances between the points ( |Rath et 
|al. ||2008| [Siitterlin et al. |2009 ) 



Besides the countless applications in time series analysis the 
use of scaling indices has been extended to the field of image 
processing for texture discrimination ( Rath &; Morfill [1 997) 



2001 Rath et al. 



and feature extraction (Jamitzky et al. 
2008) tasks. Following further this line we performed sev- 
eral non-Gaussianity studies of the CMB based on WMAP 



data using scaling indices in recent years (|Rath &; Schuecker 



2003||Rath et al. |2007||2009| |Rossmanith et al. |2009 ) 

Let us review the formalism for calculating this test statistic 
for assessing HOCs: 

In general, the SIM is a mapping that calculates for every 
point pi,i = 1, . . . , Npi X of a point set P a single value, which 
depends on the spatial position of pi relative to the group of 
other nearby points, in which the point under consideration 
is embedded in. Before we go into the details of assessing the 
local scaling properties, let us first of all outline the steps 
of generating a point set P out of observational CMB-data. 
To be able to apply the SIM on the spherical CMB data, 
we have to transform the pixelised sky S with its pixels 
at positions (0$, 0$), i = 1, N P i X , on the unit sphere to a 
point-distribution in an artificial embedding space. One way 
to achieve this is by transforming each temperature value 
T(6i,(j>i) to a radial jitter around a sphere of radius R at 
the position of the pixel center (0$, (pi). Formally, the three- 
dimensional position vector of the point pi reads as 



Xi - 

y% = 

Zi = 

with 
dR = 



(R + dR) cos(0i)sm(0i) 
(R + dR) sin(0i)sin(0i) 
(R + dR) cos(0<) 



T(0i,&)-(T) 



(i) 

(2) 
(3) 



(4) 



Hereby, R denotes the radius of the sphere while a describes 
an adjustment parameter. The mean temperature and its 
standard deviation are characterised by (T) and <tt, respec- 
tively. By the use of the normalisation we obtain for dR zero 
mean and a standard deviation of a. Both R and a should be 
chosen properly to ensure a high sensitivity of the SIM with 
respect to the temperature fluctuations at a certain spatial 
scale. For the analysis of WMAP-like CMB data, it turned 
out that this requirement is provided using R — 2 for the ra- 
dius of the sphere and coupling the adjustment parameter a 
to the value of the below introduced scaling range parameter 
r via a — r ( [Rath et al.~||2007) . Now that we obtained our 
point set P, we can apply the SIM. For every point pi we 
calculate the local weighted cumulative point distribution 
which is defined as 



^ S r (d(pi,Pj)) 



(5) 



3 = 1 



with r describing the scaling range, while s r («) and de- 
note a shaping function and a distance measure, respectively. 



The scaling index a(pi, r) is then defined as the logarithmic 
derivative of p(pi,r) with respect to r: 



*(Pi,r) 



d log p(pi, 
dlogr 



(6) 



As mentioned above, s r («) and can in general be cho- 
sen arbitrarily. For our analysis we use a quadratic gaussian 
shaping function s r (x) 
ian norm d(pi,pj) — \\j 
this specific choice of s r («) and d(*) we obtain the following 
analytic formula for the scaling indices 



i and an isotropic euclid- 

Pj\\ as distance measure. With 



£f-2(^)e-(^ 



(7) 



where we used the abbreviation dij :— d(pi,pj). As becomes 
obvious from equation |7|), the calculation of scaling indices 
depends on the scale parameter r. Therefore, we can inves- 
tigate the structural configuration in the underlying CMB- 
map in a scale-dependent manner. For our analysis, we use 
the ten scaling range parameters rk = 0.025, 0.05, 0.25, 
k = 1,2, ...10, which (roughly) correspond to sensitive l- 
ranges from Al [83; 38 7], Al = [41; 193], . . . , Al = [8; 39] 
(|Rossmanith et al. || 2009| ). 



In order to quantify the degree of agreement between the 
surrogates of different orders with respect to their signa- 
tures left in distribution of scaling indices, we calculate the 



i=l 

and the standard deviation 

N V \ ^ 2 

\2 \ 



( 1 Np 



r k ) - (a(r-fe))) 



(8) 



(9) 



of the scaling indices ol% derived from N p considered pixels 
for the different scaling ranges rk ■ N p becomes the number of 
all pixels N P i x for a full sky analysis. To investigate possible 
spatial variations of signatures of NG and to be able to mea- 
sure asymmetries we also consider the moments as derived 
from the pixels belonging to rotated hemispheres. In these 
cases the number N p of the pixels halves and their positions 
defined by the corresponding 0- and 0-intervals vary accord- 
ing to the part of the sky being considered. Furthermore, we 
combine these two test statistics by using % 2 statistics. There 
is an ongoing discussion, whether a diagonal x 2 statistic or 
the ordinary % 2 statistic, which takes into account corre- 
lations among the different random variables through the 
covariance matrix is the better suited measure. On the one 
hand it is of course important to take into account corre- 
lations among the test statistics, on the other hand it has 
been argued ( Eriksen et al.|2004| that the calculation of the 
inverse covariance matrix may become numerically unstable 
when the correlations among the variables are strong mak- 
ing the ordinary % 2 statistic sensitive to fluctuations rather 
than to absolute deviations. Being aware of this we calcu- 
lated both x 2 statistics, namely the scale dependent diago- 
nal x 2 combining the mean and the standard deviation at a 
given scale r&, and the scale-independent x 2 combining the 
mean or/and the standard deviation calculated at all scales 
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rfc, k = 1, . . . , 10 (see Rossmanith et al. (2009)). 
Further, we calculate the corresponding ordinary % 2 statis- 
tics, which is obtained by summing over the full inverse cor- 
relations matrix C _1 . In general, this is expressed by the 
bilinear form 



Full Sky 



X 2 = (M-(M)) T C - 



(10) 



where the test statistics to be combined are comprised 
in the vector M and C is obtained by cross correlating 
the elements of M. Specifically, for obtaining the scale 



dependent x full>(a 



(a(rfc)>,o- a ( rfc ) 



combining the mean and the 



standard deviation at a given scale rk the vector M 
becomes M T (Mi,M 2 ) with Mi (a(r fc )), M 2 cr a(rk) . 
Similarly, the full scale-independent x 2 statis- 
tics X f 2 uii,(c) 5 XfW« and X f 2 uii,<a), CTcK are de- 
rived from the vectors M T consisting of M T — 
(<a(ri)),...,(a(rio))), M T = (cr a(ri) , . . . , cr a(rio) ) and 
M T = ((a(ri)),...,(a(rio)),o- a ( ri ),...,o- a ( rio )), respec- 
tively. For all our investigations we calculated both % 2 
statistics and found out that the results are only marginally 
dependent from the chosen x 2 statistics. Thus, in the 
following we will only list explicit numbers for the full \ 2 
statistics, if not stated otherwise, because this measure 
yielded overall slightly more conservative results. 




Upper Hemisphere 




5 RESULTS 

To test for NGs and asymmetries in the ILC7 map and the 
NILC5 map, we compare the different surrogate maps in the 
following way: 

For each scale we calculate the mean (a(rk)) and standard 
deviation o" a ( rfc ) of the map of scaling indices a(0, </>;rfc) of 
the full sky and a set of 768 rotated hemispheres. The north- 
ern pole of the different hemispheres is located at every pixel 
centre of the full sky with N S id e — 8 in the HEALphj^] 
(Gorski et al. 2005) pixelisation scheme. The differences 



of the two classes of surrogates are quantified by the a- 
normalized deviation S 



S(Y) 



(^surro2) 



(11) 



with, Y = (tt(r/c)), cr a ( rfc ), x 2 . Every hemisphere of the set 
of 768 hemispheres delivers one deviation value S, which is 
then plotted on a sky map at that pixel position where the 
z-axis of the rotated hemisphere pierces the sky. Fig.[l]shows 
the deviations S for the mean value S({a(rk))),k = 2, 6, 10 
for the ILC7 map as derived from the comparison of the 
different classes of surrogates for the scale-independent sur- 
rogate test and for the four selected /-ranges. The following 
striking features become immediately obvious: 
First, various deviations representing features of non- 
Gaussianity and asymmetries can be found in the S'-maps for 
the ILC7 map. These features can nearly exactly be repro- 
duced when the NILC5 map is taken as input map (results 
not shown). 

Second, we find for the scale-independent surrogate test 
(leftmost column in figs. [T]) large isotropic deviations for 
the scaling indices calculated for the smallest scale shown in 



http:/ /healpix.jpl. nasa.gov/ 



Figure 2. Probability density P(a) of the first and second order 
surrogates for the scaling indices calculated for the largest scal- 
ing range no and for the /—interval Al = [2,20]. Yellow (green) 
curves denote the densities for 20 realizations of second order 
surrogates derived from the ILC7 (NILC5) map. The black lines 
are the corresponding first order surrogates. The reference frame 
for defining the upper and lower hemispheres is chosen such that 



the difference AS ■ 



Slow becomes maximal for (a) of the 



respective map and respective scale r. 



the figure. The negative values for S indicate that the mean 
of the scaling indices for the first order surrogate is smaller 
than for the second order surrogate maps. This systematic 
trend can be interpreted such that there's more structure de- 
tected in the first order surrogate than in the second order 
surrogate maps. Obviously, the random shuffle of all phases 
has destroyed a significant amount of structural information 
at small scales in the maps. 

Third, for the scale-dependent analysis we obtain for the 
largest scales (Al — [2,20]) highly significant signatures for 
non-Gaussianities and ecliptic hemispherical asymmetries at 
the largest r— values (second column in figs.[T]). These results 
are perfectly consistent with those obtained for the WMAP 5 
yr ILC map and the foreground removed maps generated by 
Tegmark, de Oliveira-Costa & Hamilton | (|2 003| o n the basis 

2009])). The only 



of the WMAP 3 year data (see |Rath et al. 
difference between this study and our previous one is that we 
now obtain higher absolute values for S ranging now from 
-4.00 < S < 3.72 for the ILC7 map and -4.36 < S < 4.50 
for the NILC5 map as compared to —3.87 < S < 3.51 for 
the WMAP 5 yr ILC map. Thus, the cleaner the map be- 
comes due to better signal-to- noise ratio and/or improved 
map making techniques the higher the significances of the 
detected anomalies, which suggests that the signal is of in- 
trinsic CMB origin. 
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Full Sky 




Upper Hemisphere 




Figure 3. Same as figure[2]but the second smallest scaling range 
r 2 and the /-interval Al = [120,300]. 



Fourth, we also find for the smallest considered scales (Al = 
[150,300]) large isotropic deviations for the scaling indices 
calculated for a small scaling range r very similar to those 
observed for the scale-independent test. 
Fifth, we do not observe very significant anomalies for the 
two other bands (Al = [20,60] and Al = [60,120]) be- 
ing considered in this study. Thus, the results obtained for 
the scale independent surrogate test can clearly be inter- 
preted as a superposition of the signals identified in the 
two /-bands covering the largest (Al — [2,20]) and small- 
est Al = [120,300]) scales. Let us investigate the observed 
anomalies in more details. We begin with a closer look at 
the most significant deviations. Fig. [2] shows the probability 
densities derived for the full sky and for (rotated) hemi- 
spheres for the scaling indices at the largest scaling range 
rio for the first and second order surrogates for the /-interval 
Al = [2, 20]. We recognize the systematic shift of the whole 
density distribution towards higher values for the upper 
hemisphere and to lower values for the lower hemisphere. As 
these two effects cancel each other for the full sky, we do no 
longer see significant differences in the probability densities 
in this case. Since the densities as a whole are shifted, the 
significant differences between first and second order surro- 
gates found for the moments cannot be attributed to some 
salient localizable features leading to an excess (e.g. second 
peak) at very low or high values in otherwise very similar 
P(a)-densities. Rather, the shift to higher (lower) values 
for the upper (lower) hemisphere must be interpreted as a 
global trend indicating that the first order surrogate map 
has less (more) structure than the respective set of second 
order surrogates. The seemingly counterintuitive result for 



Al 


Full Skv 


"Upper 
Hemisphere 


UUVV C-l 

Hemisphere 


(a(r 2 )): 


(S/%) 


(S/%) 


(S/%) 


[2, 1024] 


7.73 / > 99.8 


4.53 / >99.8 


1.87 / 96.0 


[2, 20] 
[20, 60] 
[60, 120] 
[120, 300] 


0.14 / 56.6 
0.88 / 80.6 
0.26 / 60.4 
6.97 / >99.8 


3.54 / >99.8 
1.84 / 96.4 
0.32 / 64.8 

5.36 / >99.8 


3.44 / >99.8 
1.08 / 85.2 
0.64 / 71.6 
0.92 / 83.0 










[2, 1024] 


4.16 / >99.8 


3.77 / >99.8 


0.25 / 61.8 


[2, 20] 
[20, 60] 
[60, 120] 
[120, 300] 


0.48 / 69.2 
1.70 / 95.2 
0.88 / 80.0 
3.54 / >99.8 


0.48 / 69.8 
3.18 / >99.8 
2.35 / 98.8 
1.03 / 83.4 


0.19 / 58.0 
1.02 / 84.8 
1.25 / 88.2 
3.69 / >99.8 


y 2 


(r 2 )" 






[2, 1024] 


24.55/>99.8 


14.44 / >99.8 


0.94 / 84.4 


[2, 20] 
[20, 60] 


0.90 / 85.2 
0.82 / 83.4 


7.67 / >99.8 
4.03 / 99.2 


8.47 / 99.8 
0.31 / 50.4 



[60, 120] 
[120, 300] 



0.51 / 61.4 
19.62 / >99.£ 



3.63 / 98.6 
17.17 / >99i 



1.00 / 85.2 
4.15 / 99.2 



Table 1. Deviations S and empirical probabilities p of the mean, 
standard deviation and their x 2 - com bination as derived for the 
scaling indices at the second smallest scale r 2 . The results of the 
ILC7 map are shown for the different /-bands as well as for the 
full sky and the upper and lower hemispheres. Corresponding to 
the small scale r 2 the largest values for S are calculated for small 
scale non-Gaussianities in the /-range [120, 300] and for the scale- 
independent NGs, where the phases of all Z's (Al = [2, 1024]) are 
included. 



the upper hemisphere is on the other hand consistent with a 
linear hemispherical structure analysis by means of a power 
spectrum analysis, where also a lack of power in the northern 
hemisphere and thus a pronounced hemispherical asymme- 
try was detected ( |Hansen, Banday, &; Go rski 2004] |Hansen 
et al. | |2009| . However, it has to be emphazised that the ef- 

by construction 



fects contained in the power spectrum are 
- exactly preserved in both classes of surrogates, so that the 
scaling indices measure effects that can solely be induced by 
HOCs thus being of a new, namely non-Gaussian, nature. 
Interestingly though, the linear and nonlinear hemispherical 
asymmetries seem to be correlated with each other. 
Fig. [3] is very similar to fig. [2] and shows the probability den- 
sities for the scaling indices calculated for the second small- 
est scaling range r 2 for the first and second order surrogates 
for the /-interval Al — [120,300]. The systematic shift to- 
wards smaller values for the first order surrogate for both 
hemispheres and thus for the full sky is visible. It is inter- 
esting to note that all densities derived from the ILC7 and 
NILC5 map differ significantly from each other. These dif- 
ferences can be attributed to e.g. the smoothing of the ILC7 
map. However, the systematic differences between first and 
second order surrogates induced by the phase manipulations 
prevailed in all cases - irrespective of the input map. 
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Al 


Full Skv 


"Upper 
Hemisphere 


T (fYwrpy 

-LJ u VVCl 

Hemisphere 


(a(rio)): 


(S/%) 


(S/%) 


(S/%) 


[2, 1024] 


3.75 / >99.8 


3.53 / >99.8 


1.72 / 95.4 


[2, 20] 
[20, 60] 
[60, 120] 
[120,300] 


0.64 / 74.2 
0.67/ 74.2 
0.01 / 50.5 
2.45 / 99.4 


3.24 / >99.8 
1.41 / 91.6 
2.28 / 99.0 

3.58 / >99.8 


3.41 / >99.8 
2.04 / 98.0 
2.19 / 98.6 
1.38 / 92.2 


°a(n.o) : 








[2, 1024] 


0.66 / 74.4 


3.60 / >99.8 


2.90 / >99.8 


[2, 20] 
[20, 60] 
[60, 120] 
[120,300] 


0.84 / 80.0 
2.27 / 98.6 
0.77 / 79.0 
0.60 / 73.6 


3.09 / >99.8 
2.94 / 99.8 
1.63 / 94.6 
1.61 / 95.8 


1.79 / 96.4 
0.13 / 55.0 
0.47 / 67.6 
0.81 / 79.6 


Y 2 

A («( r 10))> cr c 


*Oi )' 






[2, 1024] 


1.46 / 90.4 


9.83 / >99.8 


3.15 / 98.0 


[2, 20] 
[20, 60] 
[60, 120] 
[120,300] 


0.21 / 54.8 
2.74 / 97.2 
0.38 / 50.2 
0.26 / 57.2 


7.10 / >99.8 
5.27 / 99.6 
2.09 / 94.2 
2.23 / 96.2 


6.77 / 99.8 
0.29 / 73.6 
0.43 / 75.8 
0.19 / 60.4 



Table 2. Same as table^but for the scaling indices at the largest 
scale rio- The largest values for S are found for large scales non- 
Gaussianities in the Grange [2,20]. 



The results for the deviations \S(r)\ for the full sky and 
rotated upper and lower hemisphere are shown for all con- 
sidered /-ranges and all scales r in figs. [4]. The corresponding 
values for r<2 and no are listed in the tables [I] and [2] In table 
[3]we further summarize the results for the scale-independent 
X 2 -measures x< a >> xl a and X(a),* a - 

The main results which were already briefly discussed 
on the basis of figs. 1 become much more apparent when 
interpreting fig. [2] and tables [l] to [3] We find stable 3.7- 12a 
deviations for all r- values for S({a(rk))) and the scale- 
independent surrogate test when considering the full sky. 
This yields lvalues of S({a(r 2 ))) = 7.73 (ILC7 map) and 
S({a(r2))) = 11.06 (NILC5 map) for the scaling indices cal- 
culated for the small value r<i and S((a(rio))) = 3.75 (ILC7 
map) and S((a(rio))) = 5.77 (NILC5 map) for the largest 
radius rio- This stable r-independent effect leads to very 
high values of the deviations S for the scale-independent x 2 - 
statistics S(xf a )), where we find S(xt a \) = 5.73 (ILC7 map) 
and S(xf a )) = 27.93 (NILC5 map). It is interesting to com- 
pare these results with those obtained for the diagonal x 2 - 
statistics. In this case we find S(x 2 a )) = 57.32 (ILC7 map) 
and S(x 2 {a) ) = 119.16 (NILC5 map), which is up to an order 
of magnitude larger than the values for the full x 2 -statistics. 
These results are very remarkable, since they represent - to 
the best of our knowledge - by far the most significant detec- 
tion of non-Gaussianities in the WMAP data to date. Note 
that we used here only the mean value of the distribution of 
scaling indices, which is a robust statistics not being sensi- 
tive to contributions of some spurious outliers. Further, the 



Al 


Full Sky 


Upper 
Hemisphere 


Lower 
Hemisphere 




(S/%) 


(S/%) 


(S/%) 


[2, 1024] 


5.73 / >99.8 


9.35 / >99.8 


0.33 / 55.2 


[2, 20] 
[20, 60] 
[60 1 201 
[120, 300] 


0.97 / 95.0 
1.81 / 94.2 

1 41 / QQ 

3.17 / 92.8 


4.57 / 99.6 
2.57 / 97.4 

1 ^ / QQ 

10.53 / >99.8 


4.01 / 99.2 
2.42 / 97.0 
Q1 / 8 
1.19 / 87.8 


Y 2 • 
[2, 1024] 


5.50 / >99.8 


11.50 / >99.8 


0.66 / 79.6 


[Z, ZUJ 

[20, 60] 
[60, 120] 
[120, 300] 


U.OZ / OZ.O 

2.15 / 95.8 
1.40 /98.2 
3.10 / 99.0 


4.00 / 99.8 
3.26 / 99.4 
8.90 / >99.8 


A C\A 1 QQ f\ 

2.18 / 96.4 
2.01 / 95.6 
1.90 / 95.8 


[2, 1024] 


1.89 / 94.2 


8.38 / >99.8 


3.03 / 98.8 


[2, 20] 
[20, 60] 


0.73 / 77.4 
1.60 / 92.8 


5.64 / >99.8 
3.42 / 99.2 


6.01 / 99.8 
1.49 / 91.0 



[60, 120] 0.26 / 52.4 2.15 / 96.6 0.53 / 75.6 
[120,300] 1.68 / 92.8 5.34 / 99.8 0.22 / 63.2 



Table 3. Same as table [T] but for the scale-independent x 2 - 
statistics. Also for this statistics the largest values for S are found 
for the largest Al = [2, 20] and smallest scales Al = [120, 300] and 
for the scale-independent NGs. 

scale-independent statistics X( a ) calculated for the full sky 
represents a rather unbiased statistical approach. 
The hemisperical asymmetry for NGs on large scale (Al — 
[2,20]) finds its reflection in the results of S(r). While we 
calculate significant and stable deviations S for the upper 
and lower hemispheres separately (red and blue lines) in fig. 
[4] the results for the full sky (black lines) are not significant, 
because the deviations detected in the two hemispheres are 
complementary and thus cancel each other. Therefore, we 
obtain only for the hemispheres high values for S ranging 
from S = 3.24 up to S = 7.10 (S = 4.11 up to S = 10.82) 
for the ILC7 (NILC5) map when considering the statistics 
derived from the scaling indices for the largest scales rio 
and S = 4.01 up to S = 9.76 for the scale-independent x 2 - 
statistics. 

For the smallest scales considered so far (Al — [120,300]) 
we also find significant deviations from non-Gaussianity be- 
ing much more isotropic and naturally more pronounced at 
smaller scaling ranges r < 0.15. Thus, we obtain S — 6.97 
(ILC7 map) and S = 5.30 (NILC5 map) for S((a(r 2 ))) con- 
sidering the full sky. For the scale-independent x 2 -statistics 
the most significant signatures of NGs are detected for the 
respective upper hemispheres ranging from S = 5.16 to 
S = 10.53. To test whether all these signatures are of in- 
trinsic cosmic origin or more likely due to foregrounds or 
systematics induced by e.g. asymmetric beams or map mak- 
ing, we performed the same surrogate and scaling indices 
analysis for the five additional maps described in Section 2. 
Figs. [5] and [6] show the significance maps for the two /-ranges 



Scale- dependent non-Gaussianities in the WMAP data 9 



Al = [2,1024] Al = [2,1024] 




0.05 0.10 0.15 0.20 0.25 0.05 0.10 0.15 0.20 0.25 



Figure 4. Deviations ^(r)! for the ILC7 (left) and NILC5 (right) map as a function of the scale parameter r for the full sky (black) 
and the upper (red) and lower (blue) hemisphere. The plus signs denote the results for the mean (a(rfc)), the star-signs for the standard 
deviation <7 a ( rfc ) and the boxes for the x 2 -co m t>ination of (a(rk)) and cr a ( rk ). The shaded region indicates the 3a significance interval. 



Al = [2,20] and Al = [120,300], for which we found the 
most pronounced signatures in the ILC7 and NILC5 map. 
For the large scale NGs we find essentially the same results 
for the UILC7 map. The difference map, shows some signs of 
NGs and asymmetries, especially for large r- values. A closer 
look reveals, however, that both the numerator and denomi- 
nator in the equation for S are an order of magnitude smaller 
than the values obtained for the ILC7 (NILC5) maps. Thus 
the signal of the difference map can be considered to be 
subdominant. And even if it were not subdominant, the sig- 
nal coming from the residuals would rather diminish the 
signal in the ILC map than increase its significance, be- 
cause the foreground signal is spatially anticorrelated with 
the CMB-signal. Both the asymmetric beam map and the 
simulated coadded VW-map do not show any significant sig- 
nature for NGs and asymmetries. Finally, the simulated ILC 
map does show some signs of (galactic) north-south asym- 
metries which become smaller and therefore insignificant for 
increasing r, where we find the largest signal in the CMB 
maps. 

For the small scale NGs (Al = [120,300]) we also find that 
the UILC7-map yields similar results as the ILC7 and NILC5 
map with smaller significance. Once again the asymmetric 
beam map and the simulated coadded VW-map do not show 
significant signature for NGs and asymmetries. This is not 
the case for the simulated ILC map. Here, we find highly 
significant signatures for NGs and asymmetries, which show 
some similarities with significance patterns observed in the 
ILC7 (NILC5) map. Even much more striking features are 
detected in the difference map, where we find deviations as 
high as \S\ ~ 15 forming a very peculiar pattern in the sig- 
nificance maps for all r. One of us (G.R.) named this pattern 
'Eye of Sauron', which we think is a nice and adequate asso- 



ciation. It is worth noticing that we found the same pattern 
when analyzing other difference maps, e.g. year 7 - year 1 or 
year 2 - year 1. 

To better understand, where these features may come from 
we had a closer look at the zeroth, first and second order 
surrogate maps. It became immediately obvious that for the 
difference maps the fluctuations are systematically smaller 
in the regions in the galactic plane used for the ILC-map 
making than in the rest of the sky. This effect persists in 
the first order surrogate map and is only destroyed in the 
second order surrogates. This more (less) structure in first 
order surrogate map leads to lower (higher) values for the 
scaling indices, which can qualitatively explain the observed 
patterns in the significance maps. 

A much more detailed study of these high / effects and their 
possible origins is part of our current work but is beyond 
the scope of this paper. The results for the difference map 
shown here point, however, already towards a very interest- 
ing application of the surrogate technique. It may become a 
versatile tool to define criteria of the cleanness of maps in the 
sense of e.g. absence of artificially induced (scale-dependent) 
NGs in the map of the residual signal. Such a criterion may 
then in turn be implemented in the map making procedure 
so that ILC-like maps are not only minimizing the overall 
quadratic error in the map, but also e.g. the amount of un- 
physical NGs of the foregrounds. 



6 CONCLUSIONS 

To the best of our knowledge this work represents the first 
comprehensive study of scale-dependent non-Gaussianities 
in full sky CMB data as measured with the WMAP satellite. 
By applying the method of surrogate maps, which explic- 
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Figure 5. Deviations S((a(rk))) for the three scales r^.k = 2,6,10 (from top to bottom) for Al = [2,20]. The results are shown for 
(from left to right) the UILC7 map, the difference map 7yr ILC - 6yr ILC map, the asymmetric beam map, the coadded V and W-band 
from a standard simulation and the simulated ILC-like map (for more detailed information about the different maps see text). 





Figure 6. Same as Fig. [5] but for Al = [120,300]. Note that the scale for the color coding has significantly changed for the difference 
map (second column). 



itly relies on the scale-dependent shuffling of Fourier phases 
while preserving all other properties of the map, we find 
highly significant signatures of non-Gaussianities for very 
large scales and for the /-interval covering the first peak in 
the power spectrum. In fact, our analyses yield by far the 
most significant evidence of non-Gaussianities in the CMB 
data to date. Thus, it is no longer the question whether 
there are phase correlations in the WMAP data. It is rather 
to be figured out what the origin of these scale-dependent 
non- Gaussian signatures is. The checks on systematics we 
performed so far revealed that no clear candidate can be 
found to explain the low-/ signal, which we take to be cos- 
mological at high significance. These findings would strongly 
disagree with predictions of isotropic cosmologies with sin- 
gle field slow roll inflation. 

The picture is not that clear for the signatures found at 
smaller scales, i.e. at higher Vs. In this case we found that 
NGs can also easily be induced by the ILC map making 
procedure so that it is difficult to disentangle possible in- 
trinsic anomalies from effects induced by the preprocessing 
of the data. More tests are required to further pin down the 
origin of the detected high / anomalies and to probably un- 
cover yet unknown systematics being responsible for the low 



/ anomalies. Another way of ruling out effects of unknown 
systematics is to perform an independent observation prefer- 
ably via a different instrument as we are now able to do with 
the Planck satellite. 

In any case our study has shown that the method of sur- 
rogates in conjunction with sensitive higher order statistics 
offers the potential to become an important tool not only for 
the detection of scale-dependent non-Gaussianity but also 
for the assessment of possibly induced artefacts leading to 
NGs in the residual map which in turn may have important 
consequences for the map making procedures. 
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